classify-raster-data

Learning how to classify raster data

# load libraries
library(raster)

## Loading required package: sp

library(rgdal)

## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-3

Import CHM (Canopy Height Model)

# import TEAK CHM
chm <- raster("../NEONdata/D17-California/TEAK/2013/lidar/TEAK_lidarCHM.tif")

# test plot
plot(chm, main = "Plot TEAK CHM using raster package")

# test stretch image
image(log(chm), main = "Plot TEAK CHM image \nJust render pixels that stretch to fill the space")

Deal with 0 values

# plot hist to identify 0 values
hist(chm)

## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 32% of the raster cells were used. 100000 values used.

# set 0 values to NA
chm[chm == 0] <- NA

# look at hist again
hist(chm, xlab = "Tree height (m)")  # lots of tall trees

Read in aspect layer

# read in aspect
aspect <- raster("../NEONdata/D17-California/TEAK/2013/lidar/TEAK_lidarAspect.tif")

# plot aspect
plot(aspect, main = "Plot of TEAK aspect data")

Let’s look at how aspect relates to tree height

Are trees taller on north or south aspects?

First, we will set north-facing slopes to 1 and south-facing slopes to 2. * North-facing slopes are 0-45 and 315-360 degrees, Class = 1 * South-facing slopes are 135-225 degrees, Class = 2

Create classification matrix

# create array of values
class.m <- c(0, 45, 1,
             45, 135, NA,
             135, 225, 2,
             225, 315, NA,
             315, 360, 1)

# shape values into matrix
rcl.m <- matrix(class.m,
                ncol = 3,
                byrow = TRUE)

# let's look at our matrix
rcl.m

##      [,1] [,2] [,3]
## [1,]    0   45    1
## [2,]   45  135   NA
## [3,]  135  225    2
## [4,]  225  315   NA
## [5,]  315  360    1

Now reclassify the raster

# reclassify raster using matrix
asp.ns <- reclassify(aspect,
                     rcl.m)

# plot reclassified raster
plot(asp.ns, main = "North- and south-facing slopes")

Export Geotiff

Save reclassified matrix as output.

writeRaster(asp.ns,
            file = "../outputs/TEAK/TEAK_nsAspect2.tif",
            options = "COMPRESS=LZW", 
            NAflag = -9999)

# compress to keep file size as small as possible

Create a mask

We want to create a raster mask that will compare two rasters only at the selected points of north- and south-facing slopes. This masks out the data we are not interested in.

Important: When using a mask to compare two different rasters, need to pay attention to: Are the two rasters the same resolution? Are the locations matched up or is there georeference data in the raster?

# create a mask from asp.ns
asp.ns

## class       : RasterLayer 
## dimensions  : 577, 543, 313311  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 325963, 326506, 4102905, 4103482  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0, 2  (min, max)

# import NDVI layer
ndvi <- raster("../NEONdata/D17-California/TEAK/2013/spectrometer/veg_index/TEAK_NDVI.tif")

# plot it
plot(ndvi,
     main = "NDVI for TEAK site")

# mask data
nFacing.ndvi <- mask(ndvi, asp.ns)

# plot it
plot(nFacing.ndvi, 
     main = "NDVI for north- and south-facing slopes \nTEAK field site")